deSolve: Solving Initial Value Differential Equations in Rggplot2 and plotlyThis notebook is a demonstration of the R package deSolve. It is a direct implementation of the code in this vignette provided by the package authors, (Soetaert, Petzoldt and Setzer). However the plots are implemented using ggplot2, with plotly used for some 3-d examples. The plotting functions included in the package use base R. The purpose of this notebook is to demonstrate the use of modern plotting libraries with the output of the deSolve package.
library(deSolve)
library(tidyverse)
theme_set(theme_bw())
library(patchwork)
library(plotly)
library(metR)
We use the Lorenz model to demonstrate how to implement and solve differential equations in R. The Lorenz model describes the dynamics of three state variables, \(X\), \(Y\) and \(Z\). The model equations are: \[ \frac{dX}{dt} = a \cdot X + Y \cdot Z \] \[ \frac{dY}{dt} = b \cdot (Y - Z) \] \[ \frac{dZ}{dt} = -X \cdot Y + c \cdot Y - Z \] with the initial conditions: \[ X(0) = Y(0) = Z(0) = 1 \] where \(a\), \(b\) and \(c\) are three parameters with values of -8/3, -10 and 28 respectively.
Implementation of an IVP ODE in R can be separated in two parts: the model specification and the model application. Model specification consists of:
The model application consists of:
R functions from deSolve),R code for the Lorenz model.parameters <- c(a = -8/3, b = -10, c = 28)
state <- c(X = 1, Y = 1, Z = 1)
Lorenz <- function(t, state, parameters) {
with(as.list(c(state, parameters)),{
# rate of change
dX <- a*X + Y*Z
dY <- b * (Y-Z)
dZ <- -X*Y + c*Y - Z
# return the rate of change
list(c(dX, dY, dZ))
})
}
times <- seq(0, 100, by = 0.01)
out <- ode(y = state,
times = times,
func = Lorenz,
parms = parameters
)
head(out)
time X Y Z
[1,] 0.00 1.0000000 1.000000 1.000000
[2,] 0.01 0.9848912 1.012567 1.259918
[3,] 0.02 0.9731148 1.048823 1.523999
[4,] 0.03 0.9651593 1.107207 1.798314
[5,] 0.04 0.9617377 1.186866 2.088545
[6,] 0.05 0.9638068 1.287555 2.400161
Plot the results using ggplot2:
# Convert data to long format:
out_long <- data.frame(out) %>%
pivot_longer(cols = X:Z, names_to = "var", values_to = "value")
# Plot X, Y and Z versus time:
P1 <- out_long %>%
split(.$var) %>%
map(function(df) {
ggplot(df) +
geom_line(aes(time, value), color = "steelblue") +
scale_x_continuous(breaks = seq(0, 100, 20)) +
labs(x = "time", y = df$var[1],
title = paste0("<span style=\"font-family:'Lucida Console', monospace\">", df$var[1], "</span> versus <span style=\"font-family:'Lucida Console', monospace\">time</span>")) + # Ref: https://wilkelab.org/ggtext/
theme(panel.grid = element_blank(),
plot.title = ggtext::element_markdown())
})
# Plot Z versus X:
P1[[4]] <- ggplot(data.frame(out)) +
geom_point(aes(X, Z), shape = '.', color = "steelblue") +
ggtitle("Z versus X") +
theme(panel.grid = element_blank())
# Output plots as 2x2 grid:
(P1[[1]] + P1[[2]]) / (P1[[3]] + P1[[4]]) +
plot_annotation(title = "Figure 1: Solution of the ordinary differential equation")
Demonstrate the run times of a few different solvers from the deSolve package:
# List of integration routines from deSolve to include:
routines <-
list(
rk4 = rk4,
lsode = lsode,
lsoda = lsoda,
lsodes = lsodes,
daspk = daspk,
vode = vode
)
# Run each routine and display the run time as a data frame:
args <- list(y = state, times = times, func = Lorenz, parms = parameters)
times_list <- map(routines, function(f) system.time(exec(f, !!!args))) # https://stackoverflow.com/a/61132298
map(times_list, as.vector) %>%
dplyr::bind_rows() %>%
head(3) %>%
cbind(var = c("user", "system", "elapsed")) %>%
pivot_longer(cols = !contains("var"), names_to = "routine") %>%
pivot_wider(id_cols = routine, names_from = var)
Available methods:
rkMethod()
[1] "euler" "rk2" "rk4" "rk23" "rk23bs" "rk34f" "rk45f" "rk45ck" "rk45e" "rk45dp6" "rk45dp7" "rk78dp"
[13] "rk78f" "irk3r" "irk5r" "irk4hh" "irk6kb" "irk4l" "irk6l" "ode23" "ode45"
Define and use a new Runge-Kutta method, rKnew:
func <- function(t, x, parms) {
with(as.list(c(parms, x)),{
dP <- a * P - b * C * P
dC <- b * P * C - c * C
res <- c(dP, dC)
list(res)
})
}
rKnew <- rkMethod(ID = "midpoint",
varstep = FALSE,
A = c(0, 1/2),
b1 = c(0, 1),
c = c(0, 1/2),
stage = 2,
Qerr = 1
)
out <- ode(y = c(P = 2, C = 1),
times = 0:100, func,
parms = c(a = 0.1, b = 0.1, c = 0.1),
method = rKnew
)
head(out)
time P C
[1,] 0 2.000000 1.000000
[2,] 1 1.990000 1.105000
[3,] 2 1.958387 1.218598
[4,] 3 1.904734 1.338250
[5,] 4 1.830060 1.460298
[6,] 5 1.736925 1.580136
The function diagnostics prints several diagnostics of the simulation to the screen:
out1 <- rk4(state, times, Lorenz, parameters)
out2 <- lsode(state, times, Lorenz, parameters)
diagnostics(out1)
--------------------
rk return code
--------------------
return code (idid) = 0
Integration was successful.
--------------------
INTEGER values
--------------------
1 The return code : 0
2 The number of steps taken for the problem so far: 10000
3 The number of function evaluations for the problem so far: 40001
18 The order (or maximum order) of the method: 4
diagnostics(out2)
--------------------
lsode return code
--------------------
return code (idid) = 2
Integration was successful.
--------------------
INTEGER values
--------------------
1 The return code : 2
2 The number of steps taken for the problem so far: 12778
3 The number of function evaluations for the problem so far: 16633
5 The method order last used (successfully): 5
6 The order of the method to be attempted on the next step: 5
7 If return flag =-4,-5: the largest component in error vector 0
8 The length of the real work array actually required: 58
9 The length of the integer work array actually required: 23
14 The number of Jacobian evaluations and LU decompositions so far: 721
--------------------
RSTATE values
--------------------
1 The step size in t last used (successfully): 0.01
2 The step size to be attempted on the next step: 0.01
3 The current value of the independent variable which the solver has reached: 100.0072
4 Tolerance scale factor > 1.0 computed when requesting too much accuracy: 0
There is also a summary method for deSolve objects:
summary(out1)
Consider the Aphid model described in Soetaert and Herman (2009)1. It is a model where aphids (a pest insect) slowly diffuse and grow on a row of plants. The model equations are: \[ \frac{\partial N}{\partial t}=-\frac{\partial Flux}{\partial x} + g \cdot N \] and where the diffusive flux is given by: \[ Flux = -D \frac{\partial N}{\partial x} \] with boundary conditions \[ N_{x=0} = N_{x=60} = 0 \] and initial condition \[ N_x = 0 \text{ for } x \ne 30 \] \[ N_x = 1 \text{ for } x = 30 \] In the method of lines approach, the spatial domain is subdivided in a number of boxes and the equation is discretized as: \[ \frac{dN_i}{dt} = - \frac{Flux_{i,i+1} - Flux_{i-1,i}}{\Delta x_i} + g \cdot N_i \] with the flux on the interface equal to: \[ Flux_{i-1,i} = -D_{i-1,i} \cdot \frac{N_i - N_{i-1}}{\Delta x_{i-1,i}} \]
Define the model equations:
Aphid <- function(t, APHIDS, parameters) {
deltax <- c(0.5, rep(1, numboxes - 1), 0.5)
Flux <- -D * diff(c(0, APHIDS, 0)) / deltax
dAPHIDS <- -diff(Flux) / delx + APHIDS * r
# Return:
list(dAPHIDS)
}
Model parameters:
D <- 0.3 # m2/day diffusion rate
r <- 0.01 # /day net growth rate
delx <- 1 # m thickness of boxes
numboxes <- 60
# distance of boxes on plant, m, 1m intervals
Distance <- seq(from = 0.5, by = delx, length.out = numboxes)
Aphids are initially only present in two central boxes:
# Initial condistions: # ind/m2
APHIDS <- rep(0, times = numboxes)
APHIDS[30:31] <- 1
state <- c(APHIDS = APHIDS) # initialise state variables
The model is run for 200 days, producing output every day; the time elapsed in seconds to solve this 60 state-variable model is estimated with system.time:
times <- seq(0, 200, by = 1)
print(system.time(
out <- ode.1D(state, times, Aphid, parms = 0, nspec = 1, names = "Aphid")
))
user system elapsed
0.039 0.000 0.040
The matrix out consist of times (1st column) followed by the densities (next columns):
head(out[ , 1:5])
time APHIDS1 APHIDS2 APHIDS3 APHIDS4
[1,] 0 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
[2,] 1 1.667194e-55 9.555028e-52 2.555091e-48 4.943131e-45
[3,] 2 3.630860e-41 4.865105e-39 5.394287e-37 5.053775e-35
[4,] 3 2.051210e-34 9.207997e-33 3.722714e-31 1.390691e-29
[5,] 4 1.307456e-30 3.718598e-29 9.635350e-28 2.360716e-26
[6,] 5 6.839152e-28 1.465288e-26 2.860056e-25 5.334391e-24
Summary of the entire 1-D variable:
summary(out)
ggplot2First, a contour plot:
# Create a data frame of APHID numbers and distances for joining to the output:
Distance_APHIDS <- data.frame(
APHIDS_ = paste0("APHIDS", 1:length(Distance)),
Distance
)
# Pivot the output to longer format and join with the above data frame to add the Distance values:
out_long <- data.frame(out) %>%
pivot_longer(cols = starts_with("APHIDS"), names_to = "APHIDS_", values_to = "density") %>%
left_join(Distance_APHIDS, by = "APHIDS_")
# Plot:
ggplot(out_long) +
geom_raster(aes(time, Distance, fill = density),
interpolate = TRUE) +
geom_contour(aes(time, Distance, z = density), bins = 12,
color = "white", size = 0.1, alpha = 0.5) +
scale_fill_distiller(palette = "Spectral",
breaks = seq(0, 1, 0.2)) +
labs(x = "Time (days)", y = "Distance on plant (m)", fill = "Aphid density",
title = "Figure 2: Solution of the 1-dimensional aphid model",
subtitle = "Aphid density on a row of plants") +
guides(fill = guide_colourbar(barheight = unit(3.3, "in"),
ticks.colour = "black",
frame.colour = "black")) +
coord_cartesian(expand = F)
Line plots:
points <-
data.frame(
Distance = c(0, 10, 20, 30, 40, 50, 60),
density = c(0, 0.1, 0.25, 0.5, 0.25, 0.1, 0)
)
p1 <- out_long %>%
filter(time %in% seq(0, 200, by = 10)) %>%
ggplot() +
geom_line(aes(Distance, density, group = time, linetype = factor(time)),
color = "steelblue", show.legend = FALSE) +
geom_point(data = points, aes(Distance, density),
color = "indianred", size = 3) +
scale_x_continuous(breaks = seq(0, 60, 10)) +
scale_y_continuous(breaks = seq(0, 1, 0.2)) +
labs(x = "x", y = "", title = "Density at 21 evenly space times from 0 to 200") +
theme(panel.grid = element_blank(),
plot.title = element_text(size = 12))
p2 <- out_long %>%
filter(time == 100) %>%
ggplot() +
geom_line(aes(Distance, density, group = time, linetype = factor(time)),
color = "steelblue", show.legend = FALSE) +
geom_point(data = points, aes(Distance, density),
color = "indianred", size = 3) +
scale_x_continuous(breaks = seq(0, 60, 10)) +
scale_y_continuous(breaks = seq(0, 0.5, 0.1)) +
labs(x = "x", y = "", title = paste0("<span>Density at </span><span style=\"font-family:'Lucida Console', monospace\">time = 100</span>")) +
theme(panel.grid = element_blank(),
plot.title = ggtext::element_markdown(size = 12))
p1 + p2 + plot_annotation(title = "Figure 3: Solution of the Aphid model")
Consider the following simple DAE: \[ \frac{dy_1}{dt} = -y_1 + y_2 \] \[ y_1 \cdot y_2 = t \] where the first equation is a differential, the second an algebraic equation. To solve it, it is first rewritten as residual functions: \[ 0 = \frac{dy_1}{dt} + y_1 - y_2 \] \[ 0 = y_1 \cdot y_2 - t \]
In R we write:
daefun <- function(t, y, dy, parameters) {
res1 <- dy[1] + y[1] - y[2]
res2 <- y[2] * y[1] - t
# Return:
list(c(res1, res2))
}
yini <- c(1, 0)
dyini <- c(1, 0)
times <- seq(0, 10, 0.1)
system.time(out <- daspk(y = yini, dy = dyini, times = times, res = daefun, parms = 0))
user system elapsed
0.012 0.000 0.012
Plot using ggplot2:
pendulum <- function (t, Y, parms) {
with (as.list(Y),
list(c(u,
v,
-lam * x,
-lam * y - 9.8,
x^2 + y^2 -1
))
)
}
Consider the well-known pendulum equation:
\[ x' = u \] \[ y' = v \] \[ u' = - \lambda x \] \[ v' = -\lambda y - 9.8 \] \[ 0 = x^2 + y^2 - 1 \]
where the dependent variables are \(x\), \(y\), \(u\), \(v\) and \(\lambda\).
Implemented in R to be used with function radau this becomes:
pendulum <- function (t, Y, parms) {
with (as.list(Y),
list(c(u,
v,
-lam * x,
-lam * y - 9.8,
x^2 + y^2 -1
))
)
}
A consistent set of initial conditions are:
yini <- c(x = 1, y = 0, u = 0, v = 1, lam = 1)
Mass matrix:
M <- diag(nrow = 5)
M[5, 5] <- 0
M
[,1] [,2] [,3] [,4] [,5]
[1,] 1 0 0 0 0
[2,] 0 1 0 0 0
[3,] 0 0 1 0 0
[4,] 0 0 0 1 0
[5,] 0 0 0 0 0
Function radau requires that the index of each equation is specified; there are 2 equations of index 1, two of index 2, one of index 3:
index <- c(2, 2, 1)
times <- seq(from = 0, to = 10, by = 0.01)
out <- radau (y = yini, func = pendulum, parms = NULL, times = times, mass = M, nind = index)
Plot using ggplot2:
out_long <- data.frame(out) %>%
pivot_longer(cols = x:lam, names_to = "var")
# Plot X, Y and Z versus time:
P5 <- out_long %>%
split(.$var) %>%
map(function(df) {
ggplot(df) +
geom_line(aes(time, value), color = "steelblue") +
scale_x_continuous(breaks = seq(0, 100, 20)) +
labs(x = "time", y = df$var[1],
title = paste0("<span style=\"font-family:'Lucida Console', monospace\">", df$var[1], "</span> versus <span style=\"font-family:'Lucida Console', monospace\">time</span>")) + # Ref: https://wilkelab.org/ggtext/
theme(panel.grid = element_blank(),
plot.title = ggtext::element_markdown(size = 12))
})
# Plot y versus x:
P5[[6]] <- ggplot(data.frame(out)) +
geom_line(aes(x, y), color = "steelblue") +
ggtitle("y versus x") +
theme(panel.grid = element_blank(),
plot.title = element_text(size = 12))
# Output plots as 2x3 grid:
(P5[[4]] + P5[[5]] + P5[[2]]) / (P5[[3]] + P5[[1]] + P5[[6]]) +
plot_annotation(title = "Figure 5: Solution of the pendulum problem, an index 3 differential algebraic equation")
In this example, two state variables with constant decay are modeled:
eventmod <- function(t, var, parms) {
list(dvar = -0.1*var)
}
yini <- c(v1 = 1, v2 = 2)
times <- seq(0, 10, by = 0.1)
At time 1 and 9 a value is added to variable v1, at time 1 state variable v2 is multiplied with 2, while at time 5 the value of v2 is replaced with 3. These events are specified in a data.frame called eventdat:
ballode<- function(t, y, parms) {
dy1 <- y[2]
dy2 <- -9.8
# Return:
list(c(dy1, dy2))
}
The model is solved with ode:
out <- ode(func = eventmod, y = yini, times = times, parms = NULL, events = list(data = eventdat))
out_long <- data.frame(out) %>%
pivot_longer(cols = v1:v2, names_to = "var")
P6 <- out_long %>%
split(.$var) %>%
map(function(df) {
ggplot(df) +
geom_line(aes(time, value), color = "steelblue") +
scale_x_continuous(breaks = seq(0, 10, 2)) +
labs(x = "time", y = df$var[1],
title = paste0("<span style=\"font-family:'Lucida Console', monospace\">", df$var[1], "</span> versus <span style=\"font-family:'Lucida Console', monospace\">time</span>")) + # Ref: https://wilkelab.org/ggtext/
theme(panel.grid = element_blank(),
plot.title = ggtext::element_markdown(size = 12))
})
P6[[1]] + P6[[2]] +
plot_annotation(title = "Figure 6: A simple model that contains events")
This model describes the position (y1) and velocity (y2) of a bouncing ball:
ballode<- function(t, y, parms) {
dy1 <- y[2]
dy2 <- -9.8
# Return:
list(c(dy1, dy2))
}
An event is triggered when the ball hits the ground (height = 0) Then velocity (y2) is reversed and reduced by 10 percent. The root function, y[1] = 0, triggers the event:
root <- function(t, y, parms) y[1]
The event function imposes the bouncing of the ball:
event <- function(t, y, parms) {
y[1]<- 0
y[2]<- -0.9 * y[2]
return(y)
}
After specifying the initial values and times, the model is solved, here using lsode:
yini <- c(height = 0, v = 20)
times <- seq(from = 0, to = 20, by = 0.01)
out <- lsode(times = times, y = yini, func = ballode, parms = NULL,
events = list(func = event, root = TRUE),
rootfun = root)
data.frame(out) %>%
ggplot() +
geom_line(aes(time, height), color = "steelblue") +
ggtitle("Figure 7: A model, with event triggered by a root function") +
theme(panel.grid = element_blank())
We implement the lemming model, example 6 from (Shampine and Thompson 2000)2. Function lagvalue calculates the value of the state variable at t - 0.74. As long a these lag values are not known, the value 19 is assigned to the state variable. Note that the simulation starts at time = - 0.74:
# the derivative function
derivs <- function(t, y, parms) {
if (t < 0)
lag <- 19
else
lag <- lagvalue(t - 0.74)
dy <- r * y * (1 - lag/m)
list(dy, dy = dy)
}
# parameters
r <- 3.5; m <- 19
# initial values and times
yinit <- c(y = 19.001)
times <- seq(-0.74, 40, by = 0.01)
# solve the model
yout <- dede(y = yinit, times = times, func = derivs, parms = NULL, atol = 1e-10)
# plot using ggplot2
P8a <- ggplot(data.frame(yout)) +
geom_line(aes(time, y), color = "steelblue") +
scale_y_continuous(breaks = seq(0, 100, 20)) +
theme(panel.grid = element_blank())
P8b <- ggplot(data.frame(yout)) +
geom_path(aes(y, dy.y), color = "steelblue") +
scale_x_continuous(breaks = seq(0, 100, 20)) +
scale_y_continuous(breaks = seq(-200, 100, 100)) +
labs(y = "dy") +
theme(panel.grid = element_blank())
P8a + P8b + plot_annotation(title = "Figure 8: A delay differential equation model (Lemming model)")
We give here an example of a discrete time model, represented by a difference equation: the Teasel model as from Soetaert and Herman (2009, p287). The dynamics of this plant is described by 6 stages and the transition from one stage to another is in a transition matrix.
We define the stages and the transition matrix first:
Stages <- c("DS 1yr", "DS 2yr", "R small", "R medium", "R large", "F")
NumStages <- length(Stages)
# Population matrix
A <- matrix(nrow = NumStages, ncol = NumStages, byrow = TRUE,
data = c(0, 0, 0, 0, 0, 322.38,
0.966, 0, 0, 0, 0, 0,
0.013, 0.01, 0.125, 0, 0, 3.448,
0.007, 0, 0.125, 0.238, 0, 30.170,
0.008, 0, 0.038, 0.245, 0.167, 0.862,
0, 0, 0, 0.023, 0.75, 0)
)
The difference function is defined as usual, but does not return the “rate of change” but rather the new relative stage densities are returned. Thus, each time step, the updated values are divided by the summed densities:
Teasel <- function (t, y, p) {
yNew <- A %*% y
list (yNew / sum(yNew))
}
The model is solved using method "iteration":
out <- ode(func = Teasel, y = c(1, rep(0, 5) ), times = 0:50, parms = 0, method = "iteration")
Plot using ggplot2:
data.frame(out) %>%
pivot_longer(cols = X1:X6, names_to = "var") %>%
ggplot() +
geom_line(aes(time, value, color = var, linetype = var)) +
scale_y_continuous(breaks = seq(0, 1, 0.2)) +
ggsci::scale_color_d3(name = "Stage", labels = Stages) +
scale_linetype_discrete(name = "Stage", labels = Stages) +
ggtitle("Figure 9: A difference model solved with <span style=\"font-family:'Lucida Console', monospace\">method = \"iteration\"</span>",
subtitle = "Teasel stage distribution") +
theme(plot.title = ggtext::element_markdown(size = 12),
panel.grid = element_blank())
As an example of plotting multiple outputs on one plot we implement the simple combustion model, which can be found on http://www.scholarpedia.org/article/Stiff_systems: \[ y' = y^2 \cdot (1 - y) \]
The model is run with 4 different values of the initial conditions: \(y\) = 0.01, 0.02, 0.03, 0.04 and written to deSolve objects out1, out2, out3, out4:
combustion <- function (t, y, parms) list(y^2 * (1-y) )
yini <- 0.01
times <- 0 : 200
out1 <- ode(times = times, y = yini, parms = 0, func = combustion)
out2 <- ode(times = times, y = yini*2, parms = 0, func = combustion)
out3 <- ode(times = times, y = yini*3, parms = 0, func = combustion)
out4 <- ode(times = times, y = yini*4, parms = 0, func = combustion)
Plot with ggplot2:
map(list(out1, out2, out3, out4), data.frame) %>%
dplyr::bind_rows(.id = "out_") %>%
ggplot() +
geom_line(aes(time, X1, color = out_, linetype = out_)) +
labs(x = "time", y = "", color = "yini*i", linetype = "yini*i",
title = "Figure 10: Plotting 4 outputs in one figure") +
scale_y_continuous(breaks = seq(0, 1, 0.2)) +
ggsci::scale_color_d3() +
theme(legend.position = c(1, 0),
legend.justification = c(1, 0),
legend.background = element_blank(),
panel.grid = element_blank())
We demonstrate this with a Lotka-Volterra model, implemented in 1-D. The model describes a predator and its prey diffusing on a flat surface and in concentric circles. This is a 1-D model, solved in the cylindrical coordinate system. Note that it is simpler to implement this model in R package ReacTran.
We start by defining the derivative function:
lvmod <- function(time, state, parms, N, rr, ri, dr, dri) {
with (as.list(parms), {
PREY <- state[1:N]
PRED <- state[(N+1):(2*N)]
## Fluxes due to diffusion
## at internal and external boundaries: zero gradient
FluxPrey <- -Da * diff(c(PREY[1], PREY, PREY[N]))/dri
FluxPred <- -Da * diff(c(PRED[1], PRED, PRED[N]))/dri
## Biology: Lotka-Volterra model
Ingestion <- rIng * PREY * PRED
GrowthPrey <- rGrow * PREY * (1-PREY/cap)
MortPredator <- rMort * PRED
## Rate of change = Flux gradient + Biology
dPREY <- -diff(ri * FluxPrey)/rr/dr + GrowthPrey - Ingestion
dPRED <- -diff(ri * FluxPred)/rr/dr + Ingestion * assEff - MortPredator
return (list(c(dPREY, dPRED)))
})
}
Then we define the parameters, which we put in a list:
R <- 20 # total radius of surface, m
N <- 100 # 100 concentric circles
dr <- R/N # thickness of each layer
r <- seq(dr/2,by = dr,len = N) # distance of center to mid-layer
ri <- seq(0,by = dr,len = N+1) # distance to layer interface
dri <- dr # dispersion distances
parms <- c(Da = 0.05, # m2/d, dispersion coefficient
rIng = 0.2, # /day, rate of ingestion
rGrow = 1.0, # /day, growth rate of prey
rMort = 0.2 , # /day, mortality rate of pred
assEff = 0.5, # -, assimilation efficiency
cap = 10) # density, carrying capacity
After defining initial conditions, the model is solved with routine ode.1D:
state <- rep(0, 2 * N)
state[1] <- state[N + 1] <- 10
times <- seq(0, 200, by = 1) # output wanted at these time intervals
print(system.time(
out <- ode.1D(y = state, times = times, func = lvmod, parms = parms,
nspec = 2, names = c("PREY", "PRED"),
N = N, rr = r, ri = ri, dr = dr, dri = dri)
))
user system elapsed
0.333 0.008 0.340
The summary method provides summaries for both 1-dimensional state variables:
summary(out)
Extract the data for PREY and PRED:
prey <- subset(out, select = "PREY") %>%
data.frame() %>%
cbind(time = 0:200) %>%
dplyr::select(time, everything())
pred <- subset(out, select = "PRED") %>%
data.frame() %>%
cbind(time = 0:200) %>%
dplyr::select(time, everything())
Plot in 2-d using ggplot2:
# Create a data frame of vars (X1, X2, etc.) and r values for joining to the output:
var_r_prey <- data.frame(
var = paste0("X", 1:length(r)),
r
)
var_r_pred <- data.frame(
var = paste0("X", 101:(101+length(r)-1)),
r
)
# Pivot the output to longer format and join with the above data frame to add the Distance values:
prey_long <- data.frame(prey) %>%
pivot_longer(cols = starts_with("X"), names_to = "var") %>%
left_join(var_r_prey, by = "var")
pred_long <- data.frame(pred) %>%
pivot_longer(cols = starts_with("X"), names_to = "var") %>%
left_join(var_r_pred, by = "var")
# Plot:
P14a <- ggplot(prey_long) +
geom_raster(aes(time, r, fill = value),
interpolate = TRUE) +
geom_contour(aes(time, r, z = value), bins = 12,
color = "white", size = 0.1, alpha = 0.5) +
scale_fill_distiller(palette = "Spectral",
breaks = seq(0, 10, 1)) +
labs(x = "time", y = "", title = "PREY", fill = NULL) +
guides(fill = guide_colourbar(barheight = unit(2.7, "in"),
ticks.colour = "black",
frame.colour = "black")) +
coord_cartesian(expand = F) +
theme(plot.title = element_text(hjust = 0.5))
P14b <- ggplot(pred_long) +
geom_raster(aes(time, r, fill = value),
interpolate = TRUE) +
geom_contour(aes(time, r, z = value), bins = 12,
color = "white", size = 0.1, alpha = 0.5) +
scale_fill_distiller(palette = "Spectral",
breaks = seq(0, 10, 1)) +
labs(x = "time", y = "", title = "PRED", fill = NULL) +
guides(fill = guide_colourbar(barheight = unit(2.7, "in"),
ticks.colour = "black",
frame.colour = "black")) +
coord_cartesian(expand = F) +
theme(plot.title = element_text(hjust = 0.5))
P14a + P14b + plot_annotation(title = "Figure 14: 2-d plots of PREY and PRED") + plot_layout(guides = "collect")
Plot PREY in 3-d using plotly:
# Get the output into a list suitable for plotting a surface plt in plotly:
prey_3d <-
list(
x = times,
y = r,
z = t(subset(out, select = "PREY"))) # need to transpose this to match the x as times and y as r
# ColorBrewer Spectral palette:
my_palette <- rev(colorRampPalette(RColorBrewer::brewer.pal(9, "Spectral"))(40))
# Plot:
plot_ly(x = prey_3d$x, y = prey_3d$y, z = prey_3d$z, colors = my_palette) %>%
add_surface() %>%
layout(
scene = list(
camera=list(
eye = list(x=0, y=-2, z=1)
)
)
)
Plot PRED in 3-d using plotly:
pred_3d <-
list(
x = times,
y = r,
z = t(subset(out, select = "PRED"))) # need to transpose this to match the x as times and y as r
# Plot:
plot_ly(x = pred_3d$x, y = pred_3d$y, z = pred_3d$z, colors = my_palette) %>%
add_surface() %>%
layout(
scene = list(
camera=list(
eye = list(x=0, y=-2, z=1)
)
)
)
Consider the very simple 2-D model (100*100), containing just 1st order consumption, at a rate r_x2y2, where r_x2y2 depends on the position along the grid. First the derivative function is defined:
Simple2D <- function(t, Y, par) {
y <- matrix(nrow = nx, ncol = ny, data = Y) # vector to 2-D matrix
dY <- - r_x2y2 * y # consumption
return(list(dY))
}
Then the grid is created, and the consumption rate made a function of grid position (outer):
dy <- dx <- 1 # grid size
nx <- ny <- 100
x <- seq (dx/2, by = dx, len = nx)
y <- seq (dy/2, by = dy, len = ny)
# in each grid cell: consumption depending on position
r_x2y2 <- outer(x, y, FUN=function(x,y) ((x-50)^2 + (y-50)^2)*1e-4)
After defining the initial values, the model is solved using solver ode.2D. We use Runge-Kutta method ode45:
C <- matrix(nrow = nx, ncol = ny, 1)
ODE3 <- ode.2D(y = C, times = 1:100, func = Simple2D, parms = NULL, dimens = c(nx, ny), names = "C", method = "ode45")
Print a summary::
summary(ODE3)
Plot using ggplot2:
expand_grid(x, y) %>%
mutate(z = as.numeric(r_x2y2)) %>%
ggplot(aes(x, y, z = z)) +
geom_text_contour() + # requires the metR package
geom_contour() +
scale_x_continuous(breaks = seq(0, 100, 20)) +
scale_y_continuous(breaks = seq(0, 100, 20)) +
ggtitle("Figure 15: Contour plot of 2-D variables") +
coord_fixed()
Soetaert K, Herman PMJ (2009). A Practical Guide to Ecological Modelling. Using R as a Simulation Platform. Springer. ISBN 978-1-4020-8623-6.↩︎
Shampine L, Thompson S (2000). Solving Delay Differential Equations with dde23. URL http://www.runet.edu/~thompson/webddes/tutorial.pdf.↩︎